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Abstract 

I will review the finite density algorithm for lattice QCD based on finite 
chemical potential and summarize the associated difficulties. I will propose 
a canonical ensemble approach which projects out the finite baryon number 
sector from the fermion determinant. For this algorithm to work, it requires 
an efficient method for calculating the fermion determinant and a Monte Carlo 
algorithm which accommodates unbiased estimate of the probability. I shall 
report on the progress made along this direction with the Pade - Z2 estimator 
of the determinant and its implementation in the newly developed Noisy Monte 
Carlo algorithm. 
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1 Introduction 

Fermions at finite density or finite chemical potential is a subject of a wide range 
of interest. It is relevant to condensed matter physics, such as the Hubbard model 
away from half-filling. The research about nuclei and neutron stars at low and high 
nucleon density is actively pursued in nuclear physics and astrophysics. The subject 
of quark gluon plasma is important for understanding the early universe and is being 
sought for in relativistic heavy-ion collisions in the laboratories. Furthermore, spec- 
ulation about color superconducting phase has been proposed recently for quantum 
chromodynamics (QCD) at very high quark density |]. 

Although there are models, e.g. chiral models and Nambu Jona-Lasinio model 
which have been used to study QCD at finite quark density, the only way to study 
QCD at finite density and temperature reliably and systematically is via lattice gauge 
calculations. There have been extensive lattice calculations of QCD at finite temper- 
ature 0. On the contrary, the calculation at finite density is hampered by the lack 
of a viable algorithm. 

In this talk, I shall first review the difficulties associated with the finite density 
algorithm with chemical potentials in Sec. 2. I will then outline in Sec. 3 a proposal 
for a finite density algorithm in the canonical ensemble which projects out the nonzero 
baryon number sector from the fermion determinant. In Sec. 4, a newly developed 
Noisy Monte Carlo algorithm which admits unbiased estimate of the probability is 
described. Its application to the fermion determinant is outlined in Sec. 5. I will 
discuss an efficient way, the Pade-Z 2 method, to estimate the Tr log of the fermion 
matrix in Sec. 6. The recent progress on the implementation of the Kentucky Noisy 
Monte Carlo algorithm to dynamical fermions is presented in Sec. 7. Finally, a 
summary is given in Sec. 8. 

2 Finite Chemical Potential 

The usual approach to the finite density in the Euclidean path-integral formalism of 
lattice QCD is to consider the grand canonical ensemble with the partition function 



where the fermion fields with fermion matrix M has been integrated to give the 
determinant. U is the gauge link variable and S g is the gauge action. The chemical 
potential is introduced to the quark action with the e^ a factor in the time-forward 
hopping term and e _Ata in the time-backward hopping term. Here a is the lattice 
spacing. However, this causes the fermion action to be non-Hermitian, i.e. 75M75 7^ 
M. As a result, the fermion determinant detM[U] is complex and this leads to the 
infamous sign problem. 
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There are several approaches to avoid the sign problem: 



2.1 Fugacity Expansion 

It was proposed by the Glasgow group M that the sign problem can be circumvented 
based on the expansion of the grand canonical partition function in powers of the 
fugacity variable e^ T , 

B=3V 

ZaoQi/T,T,V) = £ e^ TB Z B (T,V), (2) 

B=-3V 

where Zb is the canonical partition function for the baryon sector with baryon number 
B. Z GG is calculated with reweighting of the fermion determinant 

_ detM\U^j\ 

ZgM ~ W[[/,or =0 - (3) 

Since the reweighting is based on the gauge configuration with /j, — 0, it avoids the 
sign problem. However, this does not work, except perhaps at small fi or near the 
finite temperature phase transition. We will dwell on this later in Sec. 3. This is 
caused by the 'overlap problem' Q where the important samples of configurations in 
the fi = simulation has exponentially small overlap with those relevant for the finite 
density. As a result, the onset of baryon begins at /i ~ m^/2 instead of the expected 
Mtv/3 which resembles the situation of the quenched approximation. 



2.2 Imaginary Chemical Potential 

In this approach, the chemical potential is taking an imaginary value \i = iv. The 
fermion determinant is real in this case and one can avoid the sign problem P, El 171. 
The partition function is 

Z GC (w/T, T, V) = Tr e -"/ T e^ T , (4) 

which is periodic with respect to v with a period of 2irT. Comparing with Eq. (|2|), 
one can in principle obtain canonical partition function Zb from the Fourier transform 

1 r 27TT 

Zb(T, V) = ^ f J Q duZ GC (zu/T, T, V)e- B ' T . (5) 

In this approach, one needs to integrate over the whole range of v from to 2nT 
after one obtains the Monte Carlo configurations of Zoo(i 1/ /T,T,V) at different v. 
In practice, it is proposed to calculate the following ratio in the two-dimensional 
Hubbard model 0, 

Z GC {iv/T,T,V) t ± „ J _ . detM(iu) 



Z GC (tv Q /T,T,V) J r v u 'detM{iv Q y 
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with a reference value vq. Several patches each centered around a different reference 
point uq are used to cover the range of v. This was successful for the two-dimensional 
Hubbard model with a 4 2 x 10 lattice up to B = 6 where the determinant was cal- 
culated exactly. While this works for a small lattice in the Hubbard model, it would 
not work for reasonably large lattices in QCD. This is because the direct calculation 
of the determinant is a V 3 (or V 2 for a sparse matrix) operation which is an imprac- 
ticable task for the quark matrix which is typically of the dimension 10 6 x 10 6 . Any 
stochastic estimation of the determinant will inevitably introduce systematic error. 
Furthermore, this will also suffer from the 'overlap' problem discussed above. Any 
Monte Carlo simulation at a reference point vq will have exponentially small overlap 
with those configurations important to a nonzero baryon density. 



2.3 Overlap Ensuring Multi-parameter Reweighting 

To alleviate the sign problem with the real chemical potential and the overlap problem 
due to reweighting, it is proposed || to do the reweighting in the multiple parameter 
space. The generic partition function Z GC in Eq. ([!]) is parametrized by a set of 
parameters a, such as the chemical potential //, the gauge coupling /?, the quark mass 
m q , etc. The partition function can be written to facilitate reweighting 

Z GC (a) = f VUdetMp, ao ] e -sAu,a o]{e ~s 3 [u, a ]+s 3 [u,a o] detM[U, a] 
J detM[U, a \ 

where the Monte Carlo simulation is carried out with the Oq set of parameters and 
the terms in the curly bracket are treated as observables. This is applied to study 
the end point in the T-/i phase diagram. In this case, the Monte Carlo simulation is 
carried out where the parameters in «o include \x = and (3 C which corresponds to the 
phase transition at temperature T c . The parameter set a in the reweighted measure 
include mu ^ and an adjusted (3 in the gauge action. The new (3 is determined 
from the Lee- Yang zeros so that one is following the transition line in the T-/i plane 
and the large change in the determinant ratio in the reweighting is compensated by 
the change in the gauge action to ensure reasonable overlap. This is shown to work 
to locate the transition line from \i = and T = T c down to the critical point on the 
4 4 and 6 3 x 4 lattices with staggered fermions ||. 

While the multi-parameter reweighting is successful near the transition line, it is 
not clear how to extend it beyond this region, particularly the T = case where one 
wants to keep the (3 and quark mass fixed while changing the \x. One still expects 
to face the overlap problem in the latter case. For large volumes, calculating the 
determinant ratio will be subjected to the same practical difficulty as discussed in the 



previous section ^2 
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3 Finite Baryon Density — A Canonical Ensemble 
Approach 



We would like to propose an algorithm to overcome the overlap problem at zero tem- 
perature which is based on the canonical ensemble approach. To avoid the overlap 
problem, one needs to lock in a definite nonzero baryon sector so that the exponen- 
tially large contamination from the zero-baryon sector is excluded. To see this, we 
first note that the fermion determinant is a superposition of multiple quark loops of 
all sizes and shapes. This can be easily seen from the property of the determinant 



detM = e^M = 1 + E (Tr\o g M) 



n=l n - 

Upon a hopping expansion of log M, Tr log M represents a sum of single loops with all 
sizes and shapes. The determinant is then the sum of all multiple loops. The fermion 
loops can be separated into two classes. One is those which do not go across the time 
boundary and represent virtual quark- ant iquark pairs; the other includes those which 
wraps around the time boundary which represent external quarks and antiquarks. The 
configuration with a baryon number one which entails three quark loops wrapping 
around the time boundary will have an energy Mb higher than that with zero baryon 
number. Thus, it is weighted with the probability e - M BN t a t compared with the one 
with no net baryons. We see from the above discussion that the fermion determinant 
contains a superposition of sectors of all baryon numbers, positive, negative and zero. 
At zero temperature where MsN t a t ^> 1, the zero baryon sector dominates and all 
the other baryon sectors are exponentially suppressed. It is obvious that to avoid the 
overlap problem, one needs to select a definite nonzero baryon number sector and stay 
in it throughout the Markov chain of updating configurations. To select a particular 
baryon sector from the determinant can be achieved by the following procedure 0: 
first, assign an U(l) phase factor e~ 1 ^ to the links between the time slices t and t + 1 
so that the link U/W is multiplied by e~ l ^/e 1 ^; then the particle number projection 
can be carried out through the Fourier transformation of the fermion determinant 
like in the BCS theory 

P N = — r dcPe-^ ! detM[(P] (9) 

Z7T JO 

where N is the net particle number, i.e. particle minus antiparticle. Note that all 
the virtual quark loops which do not reach the time boundary will have a net phase 
factor of unity; only those with a net N quark loops across the time boundary will 
have a phase factor e l ^ N which can contribute to the integral in Eq. (Q). Since QCD 
in canonical formulation does not break Z(3) symmetry, it is essential to take care 
that the ensemble is canonical with respect to triality. To this end, we shall consider 
the triality projection || 10| to the zero triality sector 



det Q M=- detM[(f> + k2ir/3]. (10) 

^ k=0,±l 
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Figure 1: Tr log M[0] for a 8 3 x 12 configuration with Wilson action as a function of 



This amounts to limiting the quark number N to a multiple of 3. Thus the triality 
zero sector corresponds to baryon sectors with integral baryon numbers. 

Another essential ingredient to circumvent the overlap problem is to stay in the 
chosen nonzero baryon sector so as to avoid mixing with the zero baryon sector 
with exponentially large weight. This can be achieved by preforming the baryon 
number projection as described above before the accept /reject step in the Monte 
Carlo updating of the gauge configuration. If this is not done, the accepted gauge 
configuration will be biased toward the zero baryon sector and it is very difficult to 
project out the nonzero baryon sector afterwords. This is analogous to the situation 
in the nuclear many-body theory where it is known Jl3[ that the variation after 
projection (Zeh-Rouhaninejad-Yoccoz method |L4], [15|]) is superior than the variation 



before projection (Peierls-Yoccoz method fll6||). The former gives the correct nuclear 
mass in the case of translation and yields much improved wave functions in mildly 
deformed nuclei than the latter. 

To illustrate the overlap problem, we plot in Fig.l Tr log M[<f>] for a configuration 
of the 8 3 x 12 lattice with the Wilson action with f3 = 6.0 and k = 0.150 which is 
obtained with 500 Z 2 noises. We see that the it is rather flat in indicating that 
the Fourier transform in Eq. (^) will mainly favor the zero baryon sector. On the 
other hand, at finite temperature, it is relatively easier for the quarks to be excited 
so that the zero baryon sector does not necessarily dominate other baryon sectors. 
Another way of seeing this is that the relative weighting factor e ~ M B^ta t can ^ e 0(1) 
at finite temperature. Thus, it should be easier to project out the nonzero baryon 
sector from the determinant. We plot in Fig. 2 a similarly obtained Tr log M[0] for a 
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8 2 x 20 x 4 $ = 4.9 n = 0.182 Wilson action lattice 
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Figure 2: Tr log M[0] for a 8 x 20 2 x 4 finite temperature configuration with dynamical 
fermion. 

configuration of the 8 x 20 2 x 4 lattice with j3 = 4.9 and k = 0.182. We see from the 
figure that there is quite a bit of wiggling in this case as compared to that in Fig. 1 
indicating that it is easier to project out a nonzero baryon sector through the Fourier 
transform at finite temperature. 

We should mention that while we think we can overcome the overlap problem and 
the determinant detM[(p] is real in this approach, nevertheless in view of the fact that 
the Fourier transform in Eq. (Q) involves the quark number N the canonical approach 
may still have the sign problem at the thermodynamic limit when iV and V are very 
large. However, we think it might work for small N such as 3 or 6 for one or two 
baryons in a finite V. This should be a reasonable start for practical purposes. 

While it is clear what the algorithm in the canonical approach entails, there are 
additional practical requirements for the algorithm to work. These include an un- 
biased estimation of the huge determinant in lattice QCD and, moreover, a Monte 
Carlo algorithm which accommodates the unbiased estimate of the probability. We 
shall discuss them in the following sections. 

4 A Noisy Monte Carlo Algorithm 

There are problems in physics which involve extensive quantities such as the fermion 
determinant which require V 3 steps to compute exactly. Problems of this kind with 
large volumes are not numerically applicable with the usual Monte Carlo algorithm 
which require an exact evaluation of the probability ratios in the accept /reject step. 
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To address this problem, Kennedy and Kuti | J proposed a Monte Carlo algorithm 
which admits stochastically estimated transition probabilities as long as they are 
unbiased. But there is a drawback. The probability could lie outside the interval 
between and 1 since it is estimated stochastically. This probability bound violation 
will destroy detailed balance and lead to systematic bias. To control the probability 
violation with a large noise ensemble can be costly. 

We propose a noisy Monte Carlo algorithm which avoids this difficulty with two 
Metropolis accept/reject steps. Let us consider a model with Hamiltonian H(U) 
where U collectively denotes the dynamical variables of the system. The major in- 
gredient of the new approach is to transform the noise for the stochastic estimator 
into stochastic variables. The partition function of the model can be written as 



Z = J[DU]e~ H ^ 

= J[DU\[Dt]P&)mt). (11) 



where /(£/,£) is an unbiased estimator of e H ^ u > from the stochastic variable £ and 
P^ is the probability distribution for £. 

The next step is to address the lower probability-bound violation. One first notes 
that one can write the expectation value of the observable O as 

(O) = J[DU}[DZ}Ps(0O(U) S ign(f)\f(U,0\/Z, (12) 

where sign(f) is the sign of the function /. After redefining the partition function to 
be 

Z = J[DU\[DZ]P e (Q\mZ)\, (13) 

which is semi-positive definite, the expectation of O in Eq. (O) can be rewritten as 

(O) = (0{U) sign(/))/(sign(/)>. (14) 

As we see, the sign of f(U,£) is not a part of the probability any more but a part 
in the observable. Notice that this reinterpretation is possible because the sign of 
f{U,£) is a state function which depends on the configuration of U and £. 

It is clear then, to avoid the problem of lower probability-bound violation, the 
accept /reject criterion has to be factorizable into a ratio of the new and old proba- 
bilities so that the sign of the estimated /(£/,£) can be absorbed into the observable. 
This leads us to the Metropolis accept/reject criterion which incidentally cures the 
problem of upper probability-bound violation at the same time. It turns out two 
accept/reject steps are needed in general. The first one is to propose updating of U 
via some procedure while keeping the stochastic variables £ fixed. The acceptance 
probability P a is 

Pa{Ui, £ - U 2 , = min(l, • (15) 
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The second accept/reject step involves the refreshing of the stochastic variables £ 
according to the probability distribution while keeping U fixed. The acceptance 
probability is 

PMti - U,&) = imn(l,l^^|) . (16) 

It is obvious that there is neither lower nor upper probability-bound violation in 
either of these two Metropolis accept/reject steps. Furthermore, it involves the ratios 
of separate state functions so that the sign of the stochastically estimated probability 
f{U,£) can be absorbed into the observable as in Eq. (14). 



Detailed balance can be proven to be satisfied and it is unbiased [12]. Therefore, 
this is an exact algorithm. 



5 Noisy Monte Carlo with Fermion Determinant 

One immediate application of NMC is lattice QCD with dynamical fermions. The 
action is composed of two parts - the pure gauge action S g (U) and a fermion action 
Sp{U) = —Tr In M(U). Both are functionals of the gauge link variables U. 

To find out the explicit form of f(U,£), we note that the fermion determinant can 



be calculated stochastically as a random walk process |TT 



tvit,m- ^ , / Tr In M . Tr In Ad . . . . , . 

TrlnM = l + TrmM(l + (1 + (...))) . (17) 



2 

This can be expressed in the following integral 

„\ oo 



g Tr In M 



~[[dr)iP v (r)i) / J[ dp T 

i=l J ° n=2 



[1 + 77}lnM77i(l + 9{p 2 - -)4 lnM?7 2 (l + 6{p 3 - -)4\nMr) 3 



where P v {r]i) is the probability distribution for the stochastic variable rji. It can be 
the Gaussian noise or the Z 2 noise (P v (r]i) = 5(\r]i\ — 1) in this case). The latter 
is preferred since it has the minimum variance [Q. p n is a stochastic variable with 
uniform distribution between and 1. This sequence terminates stochastically in finite 
time and only the seeds from the pseudo-random number generator need to be stored 



in practice. The function f(U,rj,p) ( £ in Eq. (|TTf) is represented by two stochastic 



variables r\ and p here) is represented by the part of the integrand between the the 



square brackets in Eq. (fjl|). One can then use the efficient Pade-Z 2 algorithm |19| to 
calculate the rjilnMrji in Eq. flT8]). We shall discuss this in the next section. 

Finally, there is a practical concern that Tr In M can be large so that it takes 
a large statistics to have a reliable estimate of e TrlnM from the series expansion in 
Eq. ([18]). In general, for the Taylor expansion e x = ^x n /n\, the series will start to 
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converge when x n /n\ > x n+l /(n + 1)!. This happens at n — x. For the case x = 100, 
this implies that one needs to have more than 100! stochastic configurations in the 



Monte Carlo integration in Eq. (18) in order to have a convergent estimate. Even 



then, the error bar will be very large. To avoid this difficulty, one can implement the 
following strategy. First, one notes that since the Metropolis accept/reject involves 
the ratio of exponentials, one can subtract a universal number x$ from the exponent 
x in the Taylor expansion without affecting the ratio. Second, one can use a specific 
form of the exponential to diminish the value of the exponent. In other words, one 
can replace e x with (e( x ~ X0 " N ) N to satisfy \x — Xq\/N < 1. The best choice for x$ is 
x, the mean of x. In this case, the variance of e x becomes e 5 l N — 1. 



6 The Pade — Z 2 Method of Estimating Determi- 
nants 

Now we shall discuss a very efficient way of estimating the fermion determinant 
stochastically fll9] . 

6.1 Pade approximation 

The starting point for the method is the Pade approximation of the logarithm func- 
tion. The Pade approximant to log(z) of order [K, K] at zq is a rational function 
N(z)/D(z) where deg N(z) = deg D(z) = K, whose value and first 2K derivatives 
agree with log z at the specified point zq. When the Pade approximant N(z)/D(z) is 
expressed in partial fractions, we obtain 




whence it follows 

K 

log det M = Tr logM « b TrI + £ b k ■ Tr(M + c k iy\ (20) 

k=l 

The Pade approximation is not limited to the real axis. As long as the function is 
in the analytic domain, i. e. away from the cut of the log, say along the negative real 
axis, the Pade approximation can be made arbitrarily accurate by going to a higher 
order [K, K] and a judicious expansion point to cover the eigenvalue domain of the 
problem. 
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6.2 Complex Z2 noise trace estimation 



Exact computation of the trace inverse for N x N matrices is very time consuming for 
matrices of size N ~ 10 6 . However, the complex Z 2 noise method has been shown to 
provide an efficient stochastic estimation of the trace [18, |20|, [21]]. In fact, it has been 
proved to be an optimal choice for the noise, producing a minimum variance \T2\. 

The complex Z 2 noise estimator can be briefly described as follows [|18|, p2 |. We 
construct L noise vectors rj l , t] 2 , ■ ■ ■ , t] L where rf = {rj{, rf 2 , 7/3, • ■ ■ , rf N } T , as follows. 
Each element rf n takes one of the four values {±1, ±2} chosen independently with 
equal probability. It follows from the statistics of rj 3 n that 



E[< nn>]=E[-Y,ri 



0. 



1 



B[<i ) h>]=%X)« = * 



(21) 



~ 3=1 ~ 3=1 

The vectors can be used to construct an unbiased estimator for the trace inverse of a 
given matrix M as follows: 

1 L N 



E[< rfM- 1 !] >} = E[ 



-1 



L 1 1 

j=l m,n=l 



N 



N 



E A C + (E M ™,n)[ 7 E«] 

n m^n j 

Tr M _1 . 



The variance of the estimator is shown to be f22 

Var[< ^M -1 ^ >] = E [| < ^M^T] > -Tr 



a 



1 1 2 



M 



1 N 1 N 

- T m- 1 (M^ 1 Y = - T \M~ l I 2 

/ j m,n\ m,nJ / j I m,n\ 



The stochastic error of the complex Z 2 noise estimate results only from the off- 
diagonal entries of the inverse matrix (the same is true for Z n noise for any n). 
However, other noises (such as Gaussian) have additional errors arising from diagonal 
entries. This is why the Z 2 noise has minimum variance. For example, it has been 
demonstrated on a 16 3 x 24 lattice with (3 = 6.0 and k = 0.148 for the Wilson action 
that the Z 2 noise standard deviation is smaller than that of the Gaussian noise by a 
factor of 1.54 pi . 

Applying the complex Z 2 estimator to the expression for the TrlogM. in Eq. (|2T)D . 
we find 

EWM + cO" 1 

k 

«7EEW + (M+ Cfe )-y 
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= TEEW t e fcj " ) (22) 
j fc=i 

where = (M + c/T) -1 ?]- 7 are the solutions of 

(M + c fc I)£ fcJ W, (23) 



Since M + c^I are shifted matrices with constant diagonal matrix elements, Eq. 
can be solved collectively for all values of c& within one iterative process by several 
algorithms, including the Quasi-Minimum Residual (QMR) P3[ , Multiple-Mass Min- 
imum Residual (M 3 R) |2|], and GMRES[p5|. We have adopted the M 3 R algorithm, 



which has been shown to be about 2 times faster than the conjugate gradient algo- 
rithm, and the overhead for the multiple is only 8% p6[. The only price to pay 



is memory: for each c^, a vector of the solution needs to be stored. Furthermore, 
one observes that c& > 0. This improves the conditioning of (M + c/T) since the 
eigenvalues of M have positive real parts. Hence, we expect faster convergence for 
column inversions for Eq. (F23[). 

In the next section, we describe a method which significantly reduces the stochastic 
error. 



6.3 Improved PZ estimation with unbiased subtraction 

In order to reduce the variance of the estimate, we introduce a suitably chosen set 
of traceless N x N matrices Q^, i.e. which satisfy J2n=i Qn n = 0, p = 1 ■ • • P. The 
expected value and variance for the modified estimator < ^(M -1 — J2p=i ^pQ^) 7 ? > 
are given by 

p 

E[< ^(M -1 - ]T \ p Q ip) ))v >] = Tr M" 1 , (24) 



p=i 



A M (A)=Var[<r / t(M- 1 -X:A p QW)r / >] = I £ \M~] n - £ ApQ£ re )| 2 , (25) 

p=l m^n p=l 

for any values of the real parameters X p . In other words, introducing the matrices 
Q( p ) into the estimator produces no bias, but may reduce the error bars if the are 
chosen judiciously. Further, X p may be varied at will to achieve a minimum variance 
estimate: this corresponds to a least-squares fit to the function t^M -1 // sampled at 
points rjj, j = 1 • • • L, using the fitting functions jl, 7/^^77} , p = 1 ■ ■ ■ P. 

We now turn to the question of choosing suitable traceless matrices Q^* 1 to use in 
the modified estimator. One possibility for the Wilson fermion matrix M = I — kD 
is suggested by the hopping parameter — k expansion of the inverse matrix, 

(M+cjr 1 



M + c k I (l + c.Xl-^D) 
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Ki . o K> 



This suggests choosing the matrices from among those matrices in the hopping 
parameter expansion which are traceless: 

(i + c k y 

q(2) = - d 2 

^ (1 + Cfc) 3 ^' 

r>(3) = — D 3 

Q (4) = -^(D 4 -™ 4 ), 

(1 + Cfc) 5 

qO) = - D 5 

Q (6) = y^^lD 6 -™ 6 ), 

(1 + Cfc) 7 
^.2r+l 

Q (2r+1) = 7TT ^ +2 p2r+1 . r = 3,4,5 ) ... . 

(1 + Cfc) 2r+2 

It may be verified that all of these matrices are traceless. In principle, one can include 
all the even powers which entails the explicit calculation of all the allowed loops in 
TrD 2r . In this manuscript we have only included Q*- 4 **, Q 1 - 6 -*, and Q( 2r+1 ). 



6.4 Computation of Tr log M 

Our numerical computations were carried out with the Wilson action on the 8 3 x 12 
(N = 73728) lattice with (3 = 5.6. We use the HMC with pseudofermions to generate 
gauge configurations. With a cold start, we obtain the fermion matrix Mi after 
the plaquette becomes stable. The trajectories are traced with r = 0.01 and 30 
molecular dynamics steps using k = 0.150. M 2 is then obtained from Mi by an 
accepted trajectory run. Hence Mi and M 2 differ by a continuum perturbation, and 
log[detMi/detM 2 ] ~ 0(1). 

We first calculate logdetMi with different orders of Pade expansion around zq = 
0.1 and zq = 1.0. We see from Table 1 that the 5th order Pade does not give the 
same answer for two different expansion points, suggesting that its accuracy is not 
sufficient for the range of eigenvalues of Mi. Whereas, the 11th order Pade gives 
the same answer within errors. Thus, we shall choose P[ll,ll](z) with Zq — 0.1 to 
perform the calculations from this point on. 

In Table 2, we give the results of improved estimations for Tr log Mi. We see that 
the variational technique described above can reduce the data fluctuations by more 
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Table 1: Unimproved and improved PZ estimates for log [detMx] with 100 complex 
Z2 noise vectors, k = 0.150. 



P[K,K](z) 


K = 


5 


7 


9 


11 


z = 0.1 


Original: 


473(10) 


774(10) 


796(10) 


798(10) 




Improved: 


487.25(62) 


788.17(62) 


810.83(62) 


812.33(62) 


z = 1.0 


Original: 


798(10) 


798(10) 


798(10) 


799(10) 




Improved: 


812.60(62) 


812.37(62) 


812.36(62) 


812.37(62) 



Table 2: Central values for improved stochastic estimation of log[det and rth- 
order improved Jackknife errors 5 r are given for different numbers of Z2 noise vectors. 
k, is 0.150 in this case. 



#z 2 


50 


100 


200 


400 


600 


800 


1000 


3000 


10000 


th 


802.98 


797.98 


810.97 


816.78 


815.89 


813.10 


816.53 


813.15 


812.81 


So 


±14.0 


±9.81 


±7.95 


±5.54 


±4.47 


±3.83 


±3.41 


±1.97 


±1.08 




807.89 


811.21 


814.13 


815.11 


814.01 


814.62 


814.97 






Si 


±4.65 


±3.28 


±2.48 


±1.84 


±1.50 


±1.29 


±1.12 








813.03 


812.50 


811.99 


812.86 


811.87 


812.89 


813.04 






s 2 


±2.46 


±1.65 


±1.34 


±1.01 


±0.83 


±0.72 


±0.64 








812.07 


812.01 


811.79 


812.44 


812.18 


812.99 


813.03 






S 3 


±1.88 


±1.31 


±1.01 


±0.74 


±0.58 


±0.51 


±0.44 








812.28 


812.52 


812.57 


812.86 


812.85 


813.25 


813.40 






S4 


±1.20 


±0.94 


±0.68 


±0.48 


±0.39 


±0.35 


±0.30 






5 th 


813.50 


813.07 


813.36 


813.70 


813.47 


813.54 


813.50 






s 5 


±0.82 


±0.62 


±0.47 


±0.34 


±0.29 


±0.25 


±0.22 






QtS 


813.54 


813.23 


813.22 


813.28 


813.20 


813.37 


813.26 






s 6 


±0.67 


±0.49 


±0.35 


±0.25 


±0.21 


±0.18 


±0.16 






'jts 


814.18 


813.74 


813.44 


813.42 


813.39 










s 7 


±0.44 


±0.36 


±0.26 


±0.19 


±0.16 










9 th 


813.77 


813.62 


813.49 


813.40 


813.43 










S9 


±0.40 


±0.30 


±0.22 


±0.16 


±0.14 










11 th 


813.54 


813.41 


813.45 


813.44 


813.43 










Sn 


±0.38 


±0.27 


±0.21 


±0.15 


±0.13 
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825 



820 



815 



Improved estimates for Tr log M_1 



~i r 



~i r 



Improved at kappa=.150 
Unimproved: 812.8 +/- 1.1 



810 



805 - 



800 u ' ' ' ' 1 1 1 

2 4 6 8 10 12 

order of subtraction 

Figure 3: The improved PZ estimate of TrlogM! with 50 noises as a function of the 
order of subtraction and compared to that of unimproved estimate with 10,000 noises. 
The dashed lines are drawn with a distance of 1 a away from the central value of the 
unimproved estimate. 

than an order of magnitude. For example, the unimproved error 5q = 5.54 in Table 
2 for 400 Z2 noises is reduced to Sn = 0.15 for the subtraction which includes up to 
the Q 11 matrix. This is 37 times smaller. Comparing the central values in the last 
row (i.e. the 11 th order improved) with that of unimproved estimate with 10,000 Z 2 
noises, we see that they are the same within errors. This verifies that the variational 
subtraction scheme that we employed does not introduce biased errors. The improved 
estimates of Tr log Mx from 50 Z 2 noises with errors 8 r from Table 2 are plotted in 
comparison with the central value of the unimproved estimate from 10,000 noises in 
Fig. 3. 



7 Implementation of the Kentucky Noisy Monte 
Carlo Algorithm 

We have recently implemented the above Noisy Monte Carlo algorithm to the simu- 
lation of lattice QCD with dynamical fermions by incorporating the full determinant 



directly [pjfl . Our algorithm uses pure gauge field updating with a shifted gauge cou- 
pling to minimize fluctuations in the trace log is the Wilson Dirac matrix. It gives the 
correct results as compared to the standard Hybrid Monte Carlo simulation. However, 
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the present simulation has a low acceptance rate due to the pure gauge update and 
results in long autocorrelations. We are in the process of working out an alternative 
updating scheme with molecular dynamics trajectory to include the feedback of the 
determinantal effects on the gauge field which should be more efficient than the pure 
gauge update. 

8 Summary 

After reviewing the finite density algorithm for QCD with the chemical potential, we 
propose a canonical approach by projecting out the definite baryon number sector 
from the fermion determinant and stay in the sector throughout the Monte Carlo up- 
dating. This should circumvent the overlap problem. In order to make the algorithm 
practical, one needs an efficient way to estimate the huge fermion determinant and a 
Monte Carlo algorithm which admits unbiased estimates of the probability without 
upper unitarity bound violations. These are achieved with the Pade-Z2 estimate of the 
determinant and the Noisy Monte Carlo algorithm. So far, we have implemented the 
Kentucky Noisy Monte Carlo algorithm to incorporate dynamical fermions in QCD 
on a relatively small lattice and medium heavy quark based on pure gauge updating. 
As a next step, we will work on a more efficient updating algorithm and project out 
the baryon sector to see if the finite density algorithm proposed here will live up to 
its promise. 
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